r 


A0-A056  074  SOUTHCASTCRN  MASSACHUSETTS  UNIV  NORTH  DARTMOUTH  DEPT  —ETC  F/S  8/11 - I 
AN  introduction  TO  SElf  MK  SIGNAL  PROCESSING.  (U)  H 


UUN  76  C H CHEN 
UNCLASSIFIED  EE-76-9 


AFOSR-76-2951 

AFOSR-TR-76-I096  NL 


V 


AD 

HOC  FILE  COPY  ADA056074 


afosr-tr-  7 8-  1096 


TECHNICAL  REPORT  SERIES  IN  INFORMATION  SCIENCES 
(Dr.  C.  H.  Chen,  Principal  Investigator) 


SOUTHEASTERN  MASSACHUSETTS 

UNIVERSITY 

ELECTRICAL  ENGINEERING  DEPARTMENT 


Approved  for  public  release ; 
distribution  unlimited. 


NORTH  DARTMOUTH,  MASS.  0274 

78  06  27 


Grant  AP06R  76-2951 
TR  No.  EE-78-*3 
June  8,  1978 


^ INTRODUCTION  TO 
SEISMIC  SIGNAL  PROCESSING 


Department  of  Electriced  Engineering 
Southeastern  MassMhusetts  University 
North  Deo-tfflouth,  Massachusetts  027U7 


Abstract 


This  paper  examines  the  fundamental  problem  areas  and  the  available  solutions 
in  seismic  signal  processing.  Topics  considered  include  seismic  signal  modeling, 
spectral  matching  and  the  ARMA  model,  parameter  estimation,  hcmomorphlc  versus 
predictive  deconvolution,  Kalman  filtering,  and  the  measurement  of  the  first 
aiTlval  time. 
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An  Introduction  to  Seismic  Signal  Processing 
C.  H.  Chen 

1 . Int roduct Ion 

In  recent  years,  computers  have  played  an  Increasingly  important  role  in 
seismic  studies  such  as  the  petroleim  exploration,  nuclear  detection,  earthquake 
research  and  marine  seismic  studies.  Computers  are  needed  to  process  large 
volumes  of  seismic  data  from  which  useful  information  must  be  extracted  accurately. 
A number  of  signal  processing  algorithms  have  been  developed  in  recent  years  many 
of  which  are  very  useful  for  seismic  data.  Althoxigh  the  processing  techniques 
vary  with  the  nature  of  seismic  data,  an  Important  problem  in  seismic  signal 
processing  is  deconvolution.  The  received  seismic  data  can  be  considered  as  the 
restilt  of  convolution  between  the  source  signal  and  the  transmission  medium  plus 
the  additive  Instrument  noise.  This  paper  will  be  concerned  mainly  with  the 
deconvolution  of  such  convolved  signal  as  well  as  other  seismic  signal  processing 
algorithms.  This  discussion  is  preceded  by  a study  of  seismic  signal  modeling 
as  a good  understanding  of  the  seismic  signal  generation  is  much  needed  for 
effective  signal  processing. 

2.  Seismic  Signal  Modeling 

Figures  1,  2,  and  3 depict  simplified  transmission  processes  of  seismic 
waves.  Figure  1 shows  6ui  Inhomogeneous  earth  excited  by  a deep  source.  The 
earth  is  bounded  by  two  homogeneoxis  infinite  half-spaces,  the  air  and  the 
basement  rock.  Here  the  earth  is  a distributed  pEurameter  system  governed  by 
peurtlal  differential  equations.  For  digital  processing,  the  originally  continuous 
velocity  profile  can  be  quantized  and,  as  a result,  the  earth  can  now  be  modeled 
as  a lumped  paramter  system.  If  the  time  of  signal  propagation  through  a layer 
is  short  compcured  with  the  duration  of  the  signed,  then  the  lumped  parameter 
assTimptlon  is  vedld.  By  choosing  the  depth  of  each  layer  to  be  very  small,  l.e. 
considering  many  layers,  we  can  satisfy  the  lumped  peurameter  conditions. 

To  simplify  the  anedysls,  we  can  assume  that  the  system  of  Fig.  1 is  linear 
and  time-invariant.  Let  a's  represent  the  constant  parameters  usociated  with  the 
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H layers.  Each  layer  causes  a unit  delay  for  the  seissiic  wave.  We  can  expect 


the  input  and  the  output  to  satisfy  the  linear  difference  equation, 

''n  * Vn-1  \ 

Taking  s-transform  on  both  sides  of  Eq.  (l)  we  obtain  the  ratio, 

1 


m 


a.  -1  ^ -2  ^ ^ -E 

1 V V +---+V 


= B(z) 


(1) 


(2) 


which  represents  the  transfer  function  of  an  all- pole  filter.  Our  physical  claim 
that  the  lumped  parameter  model  represents  a stable  system  is  equivalent  to  the 
mathematical  condition  that  the  transfer  function  contains  no  poles  outside  the 
unit  circle.  Eq\iation  (l)  also  represents  an  autoregressive  model  for  the 
digitized  seismic  signal. 

Figure  2 describes  the  transmission  of  teleseismic  waves.  A more  appropriate 

model  is  given  by  Fig.  3 which  shows  internal  primary  reflections  caused  by  a 

downgoing  unit  impulse  6^  applied  at  the  surface.  For  clarity,  ray  paths  are 

drawn  at  oblique  incidence  but  wave-motion  analysis  is  for  normal  incidence.  Let 

b be  the  response  of  a single  layer  with  respect  to  the  unit  impulse  input  6 . 
n n 

/■* 

By  linearity  property,  a delayed  impulse  . gives  rise  to  b and  C.'.  gives 

■ n-m  n— m n 

rise  to  Cbj^  where  C is  a constant.  By  superposition  principle,  the  total  impulse 


response  can  be  written  as 


“ c,b^  , + e«b„  « +...+  c„b  „ * € » b 

n 1 n-1  2 n-2  N n-H  n n 


(3) 


where  denotes  convolution  and  ^s  are  the  hypothetical  sources  of  strength 

given  by  the  reflection  coefficients  r at  various  layers.  Taking  the  z-transform 

n 


of  Eq.  (3)  we  have 
H(z)  * E(z)B(z) 


^ -2  ^ ^ -N 

*■1*  * ®2* 

1 + ♦ *2*”^  +...+ 


(4) 


I 


which  gives  the  transfer  function  of  a normal-incidence  reflection  seismogram 
and  is  the  ABMA  model  used  in  reflection  seismology. 


3- 


Slnce  the  layered  system  Is  assumed  to  be  both  linear  and  time  invariant, 

then  the  reflection  seismogram  y^  due  to  an  arbitrary  source  pulse  s^  is 

y = h *B_  = e *b  *s  = e *(b  *s  ) = e *<1)  (J 

•'nnTInnnn'nn'nn  ' 


where  u *=  b *s  is  defined  as  a composite  wavelet  consisting  of  the  reverberation 
n n n 

wavelet  b and  source  pulse  s . Thus  Eq.  (5)  describes  the  normal  incidence 
n n 

reflection  seismogram,  where  b^  represents  the  autoregressive  component  and 
s^*e^  the  moving  average  component.  The  basic  deconvolution  problem  is  to  filter 
y^  such  that  we  can  best  recover  the  reflection  coefficient  sequence  e^. 

3.  Spectral  Matching  and  the  ABMA  Model 

The  ARMA  model  as  derived  in  the  previous  section  is  the  most  general  linear 
seismic  signal  model.  Theoretically  speaking,  the  spectrum  of  any  physical  signal 
can  be  matched,  i.e.  fitted*  perfectly  by  an  autoregressive  model  with  an 
arbitrarily  high  order.  For  a typical  set  of  teleseismic  waveforms,  a good 
spectral  matching  based  on  the  autoregressive,  i.e.  all-pole,  model  has  been 
reported  [l][2].  If  the  ARMA,  i.e.  pole-zero,  model  is  iised,  a better  spectral 
matching  is  expected.  The  degree  of  spectred  matching  can  be  measured  by  the 
mean  squared  error  between  the  actued  and  the  modelled  signals.  Figure  U 
Illustrates  spectreul  matching  by  pole-zero  model  which  has  lower  mean  squared 
error.  Again  the  spectral  matching  is  good.  However,  the  problems  associated 
with  the  ARMA  model  are  quite  obvious:  (l)  The  order  of  the  model  must  be  finite 
in  practice.  The  linear  model  is  limited  in  its  capability  not  only  in  spectral 
matching  but  also  represents  only  a first  order  approximation  to  the  original 
signal.  (2)  Computationally  the  order  of  the  model  and  the  coefficients  in  both 
numerator  euid  denominator  must  be  determined.  This  is  far  from  being  a simple 
task.  In  fact  there  has  not  been  a satisfactory  solution  to  the  problem  of 
determining  the  numerator  polynooleQ.  and  the  order  of  the  model. 


k.  Parameter  Estimation  and  Computation 


In  this  section  we  first  consider  the  sJ.l-pole  linear  prediction  analysis 
We  assme  that  the  signal  s(n),  0 £n  £N  - 1,  csm  he  approximated  as  a weighted 
linear  summation  of  past  samples,  denoted  as 

s = - > a.  s . (6) 

n in-1 

where  a^^  are  the  predictor  coefficients  which  are  the  coefficients  of  the  AR  model, 
and  p is  the  order  of  the  filter.  The  method  of  least  squares  is  most  often  used 
[3l  to  estimate  a^  hy  minimizing  the  total  mean  squared  errors  where  the  error  is 

A 

defined  as  e = s - s . There  are  two  distinct  methods  for  estimating  the 
n n n 

parameters.  The  autocorrelation  method  minimizes  the  error  e^^  over  an  infinite 

duration.  Since  the  signal  is  of  finite  duration  in  practice,  the  infinite 

duration  signal  can  be  windowed  to  become  finite  duration.  The  autocorrelation 

matrix  is  a Toepliz  matrix  whose  special  properties  lead  to  efficient  Levinson 

and  Durbin  recursion  algorithms  for  estimating  a^.  The  second  method  is  covariance 

method  which  considers  a finite  duration  signal  only  so  it  minimizes  the  error 

over  a finite  interval.  The  covariance  matrix  is  symmetric  but  the  diagonal 

, terms  are  not  equal.  The  assumption  of  zero  values  for  data  outside  the  finite 

diiration  is  not  valid.  This  is  the  main  source  of  Inaccuracy  in  both  methods. 

The  autocorreAtlon  method  guarantees  filter  (model)  stability  while  the  covariance 
% 

method  does  not. 

Computationally  the  parameters  can  be  determined  without  major  computational 
load.  The  estimation  of  autocorrelation  or  covariance  from  data  points  may  take 
most  of  computation  time  when  the  number .of  data  points  tar  exceed  the  order  of 
the  filter,  as  is  often  the  case. 

'Hie  maximum  entropy  method  states  that  the  least  assumptions  should  be  made 
about  the  unobserved  data  points.  This  may  be  restated  by  saying  that  the  spectrum 
estimated  should  be  maximally  random  (maximum  uncertainty).  The  maximum  entropy 
solution  for.  parameters  should  be  the  same  6is  that  of  the  AR  model  except  for 
details  of  the  algorithm. 
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After  the  parameters  of  the  AR  model  are  determined  we  shall  then  use  the 
ARMA  model  by  incorporating  the  numerator  polynomials.  Let  the  denominator  be 
now  fixed,  we  can  determine  the  numerator  polynomial,  from  Eq.  (U)  by  examining 
the  ratio. 


B(z)  _ 
H(z)  “ 


-1  -2  -R 

e z + e_z  +. . .+  e z 
12  H 


(7) 


where  the  coefficients  can  be  determined  by  the  methods  which  are  used  for  the 
AR  model  as  the  inverse  z-transform  of  B(z)/H(z)  is  now  available.  Obviously  one 


iteration  may  not  be  eno\igh.  We  can  hold  e^'s  constant  and  adjust  a^^'s  by 


repeating  the  above  procedxore.  This  method  is  much  simpler  in  computation  than 
direct  estimation  [3]  of  the  parameters  of  the  ARMA  model.  We  have  tested  this 
method  on  short  length  artificial  data  sequence  to  verify  the  convergence  of  the 
recursive  procedure.  Convergence  is  verified  experimentally.  For  real  data  such 
as  of  length  102U  points,  it  will  be  difficult,  however,  to  determine  the 
coefficients  if  the  order  of  the  filter  is  high. 

Recently  the  lattice  structuires  have  been  developed  which  offer  a convenient 
visual  realization  of  the  Levinson  recursion.  For  applications  where  the  short- 
term spectrvm  changes  as  a function  of  time,  the  lattice  offers  a simple,  fast- 
converging adaptive  structure  that  has  given  results  superior  to  the  traditioned 
adaptive  transversal  filter  [•♦][5]. 

In  this  section  we  have  briefly  discussed  techniques  related  to  data  modeling 
by  least  sq\iares,  especially  the  estimation  of  ARMA  model  parameters.  The 
application  to  seismic  signal  processing  is  not  limited  to  spectral  estimation 
and  data  ccmpresslon.  There  are  good  physical  interpretation  of  xMrameters  emd 
related  quantities  such  as  the  reflection  coefficients.  The  parameters  are 
potentially  useful  features  for  classification  of  teleselsmic  events.  Good 
spectral  estimation  leads  to  accurate  computation  of  spectral  ratio  which  is 
another  useful  features  in  sesimic  discrimination.  Some  recent  articles  on 
spectral  analysis  in  seismic  data  are  [6] [71. 
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5.  Hcjmomorphlc  Versus  Predictive  Deconvolution 

In  determining  the  source  pulse  by  using  the  deconvolution  method,  both 
predictive  and  homomorphic  deconvolution  methods  have  been  extensivley  studied. 
They  represent  two  importsint  but  quite  different  approach  to  the  problem.  In 
Eq.  (5)  the  theoretical  reverberation  wavelet  b^  is  minimum  phase,  while  the 
source  pulse  s^^  is  not.  For  given  autocorrelation,  only  the  minimum  phase  pulse 
corresponsding  to  s^  can  be  determined.  The  problem  in  predictive  deconvolution 
is  thus  to  determine  an  all-pass  filter  to  obtain  a source  pulse  with  correct 
phase  characteristic.  To  do  so  an  assumption  about  the  phase  characteristic  of 
the  source  pulse  is  required.  For  the  homcsnorphic  deconvolution  [8]l9l,  the 
cutoff  quefrencies  that  produce  proper  piilse  estimate  actually  also  requires  an 
assumption  about  delay  (or  phase)  properties  of  the  source  pulse.  However, 
such  phase  assumption  is  less  critical  in  homomorphic  deconvolution  than  in 
predictive  deconvolution.  In  homomorphic  deconvolution  an  exponentieO.  weighting 
of  the  data  sequence  is  usually  necessary  to  remove  ccanputational  instability  due 
to  the  nonminimm  phase  source  pulse.  Figure  5 shows  some  results  of  cepstral 
analysis  on  the  teleseismic  records. 

In  the  homomorphic  deconvolution  the  reflection  coefficient  sequence  can  be 
determined  once  the  complex  cepstnan  corresponding  to  the  source  pulse  and 
reverberation  is  removed.  For  the  predictive  deconvolution,  one  assumes  that 


the  reflection  coefficients  are  a random  uncorrelated  sequence  to  be  estimated. 

6.  Kalman  Filtering  Approach 

The  Kalman  filtering  approach  can  be  used  to  obtain  optimal  smoothed  estimates 


of  the  reflection  coefficient  sequence  from  seismic  traces  with  noise  [10].  The 

seismic  trace  can  be  interi>reted  as  the  sum  of  additive  noise  and  the  output  of 

y 

a linear  system,  with  rec^nse  u given  by  Eq.  (5),  excited  by  white  noise 

n 


y 

corresponding  to  th6  reflection  coefficient  sequence.  So  the  estimation  of 


reflection  coefficient  is  now  the  same  problem  as  estimating  the  random 
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disturbance  in  a state  equation.  Assmption  is  made  that  basic  composite  wavelet 
Is  known  a priori,  and  that  the  both  additive  noise  and  the  white  noise  for  the 
reflection  coefficient  sequence  have  known  covariance  matrices,  either  of  which 
may  be  time  veurylng.  The  Kalman  filtering  approach  permits  a more  flexible 
modeling  assumption  than  the  Wiener  filtering  to  the  predictive  deconvolution. 

7.  Measuring  the  First  Arrived.  Time 

In  microearthquakes  and  in  explosive  events,  for  example,  it  will  be  necessary 
to  measure  accurately  the  first  arrived,  time.  The  ambiguity  associated  with  the 
meas\irement  of  the  first  arrived,  time  is  due  to  the  facts  that  the  signal  is 
contaminated  by  noise  and  the  wave  shape  of  the  first  arrival  is  tinknown.  Thus 
the  methods  of  "beam  forming"  emd  "matched  filtering"  eure  not  acceptable. 

However,  the  first  eurrival  does  occur  in  a context,  the  features  of  which  cem  be 
determined  more  reliably.  Anderson  [11]  developed  a simple  but  robust  algorithm 
for  automatic  anedysis  of  microearthquake  data  to  pick  the  first  arrival  with 
good  accuracy.  The  algorithm  uses  informations  frcan  multiple  levels  such  an  a 
tentative  location,  detection  of  the  event,  arrival  time  residuals,  and  the 
features  which  represent  physical  measuranents  of  the  seismic  wave  including  the 
first  and  second  zero  crossing,  the  first  maximum  in  the  half  cycle,  etc.  Syntactic 
pattern  recognition  should  be  a useful  approach  to  Improve  the  algorithm. 
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